##############################################
### R code for PS 659 replictation project
### Spring 2013 (Apr. 30th)
###############################################

#working directory
setwd("C://Program Files/R/data/")

#load libraries
library(foreign) # to read non-R dataset
library(Zelig)
library(xlsxjars) #required for xlsx package
library(rJava)   #required for xlsx package
library(xlsx) #reading excel files
library(MASS) #probit model (polr)
library(arm) # for automatic coefficient plot
library(xtable) #for producing the output as latex code
library(separationplot) #for separation plot
library(car)

### Read datasets 
#Chapman & Roeder data (2007, APSR)
chapman<-read.xlsx("PartitionPosting1.xlsx", 2, keepFormulas=FALSE, as.data.frame=TRUE)
head(chapman)
d1<-chapman
summary(d1)

#Sambanis & Schulhofer-Wolh data (2009, International Security)
sambanis<-read.dta("IS2009_replicdata.dta")
head(sambanis)
d2<-sambanis
summary(d2)

#Fearon & Laitin data (2003, APSR)
fearon<-read.dta("fearon.dta")
head(fearon)
d3<-fearon
summary(d3)


##Cammet and Malesky data -- only the negotiated cases
#cammet<-read.dta("cammet.dta")
#head(cammet)
#d4<-cammet
#summary(d4)
##clean up the dta format (has NA rows)
#d4<-d4[1:53,]



###Fortna data
fortna<-read.dta("fortna.dta")
head(fortna)
d5<-fortna
summary(d5)

#######################Exactly replicate the model in CR (model 4 of table 3)

m1<- glm(SurvivalofPeace ~ Partition + Separation + Autonomy + Pre.warDemocracy + WarDuration + WarDeaths + ArmedForces + GDPperCapita + PeaceOperations,
        data=d1, family=binomial(link=probit))

summary(m1, digits=2)
coefplot(m1, main="Survival of Peace - Binary (Chapman and Roeder, 2007)", 
        cex.var=0.5)


#########################Exactly replicate the model in SS (table 4)

m2 <- glm(norecur2 ~ part + logcost + factnum + anypko + treaty + isxp2 + ef + lnmaddpre_i + imaddgro,
          data=d2, family=binomial(link=probit))
summary(m2)
coefplot(m2, main="", cex.var=0.7)

m3<-glm(norecur2_v2 ~ part + logcost + factnum + anypko + treaty + isxp2 + ef + lnmaddpre_i + imaddgro,
  data=d2, family=binomial(link=probit))
summary(m3)
coefplot(m3, main="", cex.var=0.7)

m4<-glm(warnov2_01 ~ part + logcost + factnum + anypko + treaty + isxp2 + ef + lnmaddpre_i + imaddgro,
  data=d2, family=binomial(link=probit))
summary(m4)
coefplot(m4, main="", cex.var=0.7)

m5<-glm(norecur2 ~ part2 + logcost + factnum + anypko + treaty + isxp2 + ef + lnmaddpre_i + imaddgro,
  data=d2, family=binomial(link=probit))
summary(m5)
coefplot(m5, main="", cex.var=0.7)

m6<-glm(norecur2_v2 ~ part2 + logcost + factnum + anypko + treaty + isxp2 + ef + lnmaddpre_i + imaddgro,
  data=d2, family=binomial(link=probit))
summary(m6)
coefplot(m6, main="", cex.var=0.7)

m7<-glm(warnov2_01 ~ part2 + logcost + factnum + anypko + treaty + isxp2 + ef + lnmaddpre_i + imaddgro,
  data=d2, family=binomial(link=probit))
summary(m7)
coefplot(m7, main="", cex.var=0.7)


#Extensions

m8<- glm(SurvivalofPeace ~ Partition + Pre.warDemocracy + WarDuration + WarDeaths + ArmedForces + GDPperCapita + PeaceOperations,
         data=d1, family=binomial(link=probit))

summary(m8, digits=2)
coefplot(m8, main="Partition on Post-war Peace", 
         cex.var=0.5)



m9<- glm(SurvivalofPeace ~ Partition,
         data=d1, family=binomial(link=probit))

summary(m9, digits=2)
coefplot(m9, main="Partition on Post-war Peace", 
         cex.var=0.5)


m10<-glm(warnov2_01 ~ part2,
        data=d2, family=binomial(link=probit))   #works at 90% significance
summary(m10)
coefplot(m10, main="", cex.var=0.7)     


#Daga merge
yrst<-d1$YearBegin+1900
d1<-cbind(yrst,d1)

dd<-merge(d1,d2, by=c("clust2","yrst"), all.x=T, all.y=T)


#combine three measures of partition (Partition + Separation + Autonomy)


all3<-dd$Partition
all3[dd$Separation==1]<-1
all3[dd$Autonomy==1]<-1
table(all3)
table(dd$Partition)
table(dd$Separation)
table(dd$Autonomy)


dd<-cbind(dd,all3)

#combine two measures of partition (Partition + Separation)


partsep<-dd$Partition
partsep[dd$Separation==1]<-1
table(partsep)
dd<-cbind(dd,partsep)


#create the unitarism variable (-(Participationi+Separation))

uni<-rep(NA,length(partsep))
uni[partsep==0]<-1
uni[partsep==1]<-0
table(partsep)
table(uni)


#see how the two variables match!

sum1<-sum2<-sum3<-sum4<-sum5<-sum6<-rep(0,184)

sum1[dd$SurvivalofPeace==dd$norecur2]<-1     #comparing the depedent variable 1
sum2[dd$SurvivalofPeace==dd$norecur2_v2]<-1    #comparing the dependent variable 2
sum3[dd$SurvivalofPeace==dd$warnov2_01]<-1   #comparing the dependent variable 3 -- residual violence
sum4[dd$Partition==dd$part] <-1     #comparing the independent variable -- partition (loose)
sum5[dd$Partition==dd$part2]<-1      #comparing the independent variable -- partition (strict)
sum6[dd$partsep==dd$part2]<-1    #comparing the independent variable -- partition (partsep & strict)

table(sum1)
table(sum2)
table(sum3)
table(sum4)
table(sum5)
table(sum6)

dd<-cbind(dd,sum1,sum2,sum3,sum4,sum5,sum6)

table(dd$Partition)
table(dd$part)
table(dd$part2)



## dependent var from RS, independent var from RS

m11<- glm(SurvivalofPeace ~ partsep,
         data=dd, family=binomial(link=probit))      # Report!

summary(m11, digits=2)
coefplot(m11, main="Partition on Post-war Peace", 
         cex.var=0.5)

## dependent var from SS, independent var from RS

m12<-glm(warnov2_01 ~ partsep,
         data=dd, family=binomial(link=probit))
summary(m12)
coefplot(m6, main="", cex.var=0.7)


m13<-glm(norecur2 ~ partsep,
         data=dd, family=binomial(link=probit))

summary(m13)
coefplot(m7, main="", cex.var=0.7)


m14<-glm(norecur2_v2 ~ partsep,
         data=dd, family=binomial(link=probit))

summary(m14)
coefplot(m7, main="", cex.var=0.7)



## dependent var from RS, independent var from SS

m15<-glm(SurvivalofPeace~ part,
         data=dd, family=binomial(link=probit))
summary(m15)
coefplot(m6, main="", cex.var=0.7)


m16<-glm(SurvivalofPeace ~ part2,
         data=dd, family=binomial(link=probit))

summary(m16)
coefplot(m7, main="", cex.var=0.7)


#### Compare the partition variable
comparison <-cbind(dd$cname, name=dd$conflict, CR=dd$partsep, SS=dd$part2, sum5)


#modeling on the cases agreed upon

agreed.partition<-rep(NA,184)
agreed.partition[sum6==1]<-1
table(agreed.partition)
agreed.data<-cbind(agreed.partition, partition=dd$partsep, survivalofpeace=dd$SurvivalofPeace)
agreed.data<-agreed.data[complete.cases(agreed.data), ]
agreed.data<-as.data.frame(agreed.data)
dim(agreed.data)

m17<-glm(survivalofpeace ~ partition,
         data=agreed.data, family=binomial(link=logit))    #Report!!

summary(m17)
coefplot(m17, main="", cex.var=0.7)

#Data with cases agreed both in terms of participation and recurrence

agreed.recurrence<-rep(NA,184)
agreed.recurrence[sum1==1]<-1 
agreed.data2<-cbind(agreed.partition, agreed.recurrence, partition=dd$part2, survivalofpeace=dd$norecur2_v2)
agreed.data2<-agreed.data2[complete.cases(agreed.data2), ]
agreed.data2<-as.data.frame(agreed.data2)
dim(agreed.data2)


m18<-glm(survivalofpeace ~ partition,
         data=agreed.data2, family=binomial(link=logit))
summary(m18)
coefplot(m18, main="", cex.var=0.7, add=TRUE)


#since it's all about coding, I am cross checking the data with Page's data
clust2<-dd$clust2
clust2<-as.data.frame(clust2)
ref<-cbind(clust2, yrst=dd$yrst, sum1=dd$sum1, survivalofpeace.rc=dd$SurvivalofPeace, survivalofpeace.ss=dd$norecur2, sum6=dd$sum6, partition.rc=dd$partsep, partition.ss=dd$part2)
table(ref$sum1)


########### Conclusion: Partition is not necessarily a superior solution.

## produce graphs to summarize the results!


### Descriptive statistics

#1. bar plot
count1<-table(dd$SurvivalofPeace)
#count1<-c(41,31)
count2<-table(dd$norecur2)
#count2<-c(22,109)
count3<-table(dd$Partition)
#count3<-c(65,7)
count4<-table(dd$part2)
#count4<-c(129,22)
par(mfrow=c(2,2))



display.brewer.all()
display.brewer.all(type="div")
Sys.sleep(2)
display.brewer.all(type="seq")
Sys.sleep(2)
display.brewer.all(type="qual")

#Color scheme
par(mfrow=c(2,2))
display.brewer.pal(8,"Accent")
accent <- brewer.pal(8,"Accent")
display.brewer.pal(11,"Spectral")
spectral <-brewer.pal(11,"Spectral")
greens<-brewer.pal(6,"Greens")
display.brewer.pal(6,"Greens")
purple<-brewer.pal(6,"BuPu")
display.brewer.pal(6,"BuPu")
blues<-brewer.pal(6,"Blues")
set3<-brewer.pal(12,"Set3")
paired<-brewer.pal(12,"Paired")



barplot(count1,main="Dependent variable (C&R)",
        xlab="Recurrence of war/ Peace", ylab="Frequency", col=spectral[2])
barplot(count2,main="Dependent variable (S&S)",
        xlab="Recurrence of war/ Peace", ylab="Frequency", col=spectral[10])
barplot(count3,main="Independent variable (C&R)",
        xlab="Non-partition/ Partition", ylab="Frequency", col=spectral[2])
barplot(count4,main="Independent variable (S&S)",
        xlab="Non-partition/ Partition", ylab="Frequency", col=spectral[10])

#Matching graph
par(mfrow=c(2,1))

base<-rep(1:184)
#base<-as.vector(base)
#sum1<-as.vector(sum1)
plot(base,dd$sum1,type="h", col=accent[6], lwd=2,
     main="Matching the dependent variable",
     xlab="Peace (sorted) ", ylab="Matched", 
     axes=F)

plot(base,dd$sum6,type="h", col=accent[7], lwd=2,
     main="Matching the independent variable",
     xlab="Partition (sorted) ", ylab="Matched",
     axes=F)



# Coefficient plot (the original model)

par(mfrow=c(1,2))
coefplot(m1, col=spectral[2], main="CR original Model", cex.var=0.5)
coefplot(m5, col=spectral[10], main="SS original Model", cex.var=0.7)


## Coefficient plot 2 (extended/replicated models)

par(mfrow=c(2,1))
coefplot(m17, col=accent[6], main="(M1) Cases only with agreed partition", cex.var=0.9)
coefplot(m18, col=accent[7], main="(M2) Cases only with agreed partition and peace", cex.var=0.9)







#### 2. Predicted Probability & spagetti plot -- Model 17 (model with agreed partition)
data<-agreed.data
sophie3<-m17
set.seed(123)


lwgvec <- seq(min(data$partition), max(data$partition), length.out=100)

predvec <- NULL
for(i in 1:length(lwgvec)){
  values <- c(1, lwgvec[i])
  pred <- sophie3$coefficients %*% values
  predvec <- c(predvec, pred)
}

predprob <- 1/(1+exp(-predvec))
#predprob <- -predvec



### WITH UNCERTAINTY 

draw <- mvrnorm(100, sophie3$coefficients, vcov(sophie3))

predmat <- matrix(NA, nrow=100, ncol=100)
for(j in 1:100){
  predvec <- NULL
  for(i in 1:length(lwgvec)){
    values <- c(1, lwgvec[i])
    pred <- draw[j,] %*% values
    predvec <- c(predvec, pred)
  }
  predmat[j,] <- 1/(1+exp(-predvec))
}

# plotting confidence interval (spaghetti plot type)
# Box


par(mfrow=c(1,1))
plot(lwgvec, predprob, type="n", ylim=c(0,1), col=accent[5], axes=F,
     ylab="Survival of Peace", xlab="Partition", main="(M1) Predicted Probability")
box()
#X axis
axis(1, at=0, lab="Non-partition")
axis(1, at=1, lab="Partition")
#Y axis
axis(2, at=0.1, lab="Recurrence")
axis(2, at=1, lab="Peace")



for(i in 1:100){
  points(lwgvec, predmat[i,], type="l", lwd=0.2, col=spectral[4])
}
lines(lwgvec, predprob, type="l", col=spectral[2], lwd=2)



#Rug at the top (placing rug at the top side=3, defulat;bottom)
rug(data$partition)





#### 3. Predicted Probability & spagetti plot (model 18, agreed partition & peace)
data<-agreed.data2
sophie3<-m18
set.seed(123)


lwgvec <- seq(min(data$partition), max(data$partition), length.out=100)

predvec <- NULL
for(i in 1:length(lwgvec)){
  values <- c(1, lwgvec[i])
  pred <- sophie3$coefficients %*% values
  predvec <- c(predvec, pred)
}

predprob <- 1/(1+exp(-predvec))
#predprob <- -predvec



### WITH UNCERTAINTY 

draw <- mvrnorm(100, sophie3$coefficients, vcov(sophie3))

predmat <- matrix(NA, nrow=100, ncol=100)
for(j in 1:100){
  predvec <- NULL
  for(i in 1:length(lwgvec)){
    values <- c(1, lwgvec[i])
    pred <- draw[j,] %*% values
    predvec <- c(predvec, pred)
  }
  predmat[j,] <- 1/(1+exp(-predvec))
}

# plotting confidence interval (spaghetti plot type)
# Box


par(mfrow=c(1,1))
plot(lwgvec, predprob, type="n", ylim=c(0,1), col=accent[5], axes=F,
     ylab="Survival of Peace", xlab="Partition", main="(M2) Predicted Probability")
box()
#X axis
axis(1, at=0, lab="Non-partition")
axis(1, at=1, lab="Partition")
#Y axis
axis(2, at=0.1, lab="Recurrence")
axis(2, at=1, lab="Peace")



for(i in 1:100){
  points(lwgvec, predmat[i,], type="l", lwd=0.2, col=spectral[8])
}
lines(lwgvec, predprob, type="l", col=spectral[10], lwd=2)



#Rug at the top (placing rug at the top side=3, defaulat;bottom)



############### Finally, two scenarios!


##################################
##### probability of two scenarios   Model 17 (m17)
#################################


##### Generate the probability of war based on two different scenarios(partition/-partition)

data<-agreed.data
sophie3<-m17
# setting up vectors to be filled
outmin <- outmin <- NULL  
set.seed(123)
draw <- mvrnorm(n=1000,sophie3$coefficients,vcov(sophie3))   
dim(draw)

# the vector of values: 1 for intercept
effectmin <- c(1,0)   
effectmax <- c(1,1)    

outmin <- draw%*%effectmin
outmax <- draw%*%effectmax

dim(outmin)
dim(outmax)

# in log scale (since it's a logit model..)
outmin<- 1/(1+exp(-outmin))
dim(outmin)
summary(outmin)

outmax <- 1/(1+exp(-outmax))
dim(outmax)
summary(outmax)

plot(density(outmin))
plot(density(outmax)) 

######### Now graphing!

#Prepare the box (outside)
plot(density(outmin), type="n", ylim=c(0,5), xlim=c(0,1), axes=F,
     ylab="Density", xlab="Probability of Peace", main="(M1)Two Scenarios of Partition and -Partition on Peace")
box()
axis(1, at=0, lab=0)
axis(1, at=0.5, lab=0.5)
axis(1, at=1, lab=1)


d<-density(outmin) #
s<-density(outmax) #

lines(d, col=spectral[11], lwd=0.8) 
polygon(d, col=spectral[9], density=30) 
abline(v=mean(outmin),b=0, col=spectral[11], lwd=0.7, lty=2)

lines(s, col=spectral[1], lwd=0.8) 
polygon(s, col=spectral[3], density=30, opacity=0.5)
abline(v=mean(outmax), b=0, col=spectral[1], lwd=0.7, lty=2)

#Now legend!
legend("topleft",c("Partition","Non-partition"),col=c(spectral[3],spectral[9]),lty=c(1,1), lwd="3")
plot(s)



##################################
##### probability of two scenarios   Model 18 (m18)
#################################


##### Generate the probability of war based on two different scenarios(partition/-partition)

data<-agreed.data2
sophie3<-m18
# setting up vectors to be filled
outmin <- outmin <- NULL  
set.seed(123)
draw <- mvrnorm(n=1000,sophie3$coefficients,vcov(sophie3))   
dim(draw)

# the vector of values: 1 for intercept
effectmin <- c(1,0)   
effectmax <- c(1,1)    

outmin <- draw%*%effectmin
outmax <- draw%*%effectmax

dim(outmin)
dim(outmax)

# in log scale (since it's a logit model..)
outmin<- 1/(1+exp(-outmin))
dim(outmin)
summary(outmin)

outmax <- 1/(1+exp(-outmax))
dim(outmax)
summary(outmax)

plot(density(outmin))
plot(density(outmax)) 

######### Now graphing!

#Prepare the box (outside)
plot(density(outmin), type="n", ylim=c(0,4), xlim=c(0,1), axes=F,
     ylab="Density", xlab="Probability of Peace", main="(M2) Two Scenarios of Partition and -Partition on Peace")
box()
axis(1, at=0, lab=0)
axis(1, at=0.5, lab=0.5)
axis(1, at=1, lab=1)


d<-density(outmin) #
s<-density(outmax) #

lines(d, col=spectral[11], lwd=0.8) 
polygon(d, col=spectral[9], density=30) 
abline(v=mean(outmin),b=0, col=spectral[11], lwd=0.7, lty=2)

lines(s, col=spectral[1], lwd=0.8) 
polygon(s, col=spectral[3], density=30, opacity=0.5)
abline(v=mean(outmax), b=0, col=spectral[1], lwd=0.7, lty=2)

#Now legend!
legend("topleft",c("Partition","Non-partition"),col=c(spectral[3],spectral[9]),lty=c(1,1), lwd="3")
plot(s)





## diagnostics
outlierTest(m17)
par(mfrow=c(2,2))
plot(m17)
plot(m18)

#########################################################
########### Second approach
### cross checking with Fortna's data
#########################################################

fortna2<-fortna10<-ref$survivalofpeace.ss
fortna2[24]<-1
fortna10[19]<-0
fortna10[21]<-0
fortna10[45]<-0
fortna10[62]<-0
fortna10[77]<-0
fortna10[94]<-0
fortna10[117]<-0
table(fortna10)
table(fortna2)
table(ref$survivalofpeace.ss)


table(agreed.partition)
checked.data<-cbind(agreed.partition, partition=ref$partition.ss, 
                    survivalofpeace.ss=ref$survivalofpeace.ss, 
                    fortna2=fortna2, fortna10=fortna10)
checked.data<-checked.data[complete.cases(checked.data), ]
checked.data<-as.data.frame(checked.data)
dim(checked.data)

# descriptive statistics

table(checked.data$partition)
ch<-rep(1:18)
ch.2<-rep(1,18)
ch.10<-rep(1,18)
ch.10[1]<-0
ch.10[2]<-0
ch.10[8]<-0
ch.10[10]<-0
ch.10[13]<-0
ch.10[17]<-0
ch.10[18]<-0
ch.cr<-rep(0,18)
ch.cr[3]<-1
ch.ss<-rep(1,18)
ch.ss[3]<-0
ch.10.cr<-c(1,1,1,0,0,0,0,1,0,1,0,0,1,0,0,0,1,1)
ch.10.ss<-c(0,0,0,1,1,1,1,0,1,0,1,1,0,1,1,1,0,0)
ch.10.cr[ch.10.cr==0]<-NA
ch.10.ss[ch.10.ss==0]<-NA
ch.cr[ch.cr==0]<-NA
ch.ss[ch.ss==0]<-NA
#Matching graph
par(mfrow=c(2,1))

plot(ch,ch.cr,type="h", col=spectral[1], lwd=20,
     main="Cross checking 'peace' with Fortna data (2yr)",
     axes=F, xlim=c(1,23), ylim=c(0,1),
     xlab="", ylab="Matched")
points(ch, ch.ss, type="h", col=spectral[11], lwd=20)
legend("right",c("Matches CR data","Matches SS data"),
       col=c(spectral[1],spectral[11]),lty=c(1,1), lwd="4", cex=0.7)


plot(ch,ch.10.ss,type="h", col=spectral[1], lwd=20,
     main="Cross checking 'peace' with Fortna data (10yr)",
     axes=F, xlim=c(1,23), ylim=c(0,1),
     xlab="", ylab="Matched")
points(ch, ch.10.cr, type="h", col=spectral[11], lwd=20)
legend("right",c("Matches CR data","Matches SS data"),
       col=c(spectral[1],spectral[11]),lty=c(1,1), lwd="4", cex=0.7)



#### Run the model (M3) - cross checked peace
m19<-glm(fortna2~partition,
         data=checked.data, family=binomial(link=logit))
summary(m19)
m20<-glm(fortna10~partition,
         data=checked.data, family=binomial(link=logit))
summary(m20)


## Coefficient plot 3 (extended/replicated models)

par(mfrow=c(2,1))
coefplot(m19, col=accent[5], main="(M3) Agreed partition and checked peace (2yr)", cex.var=0.9)
coefplot(m20, col=accent[8], main="(M4) Agreed partition and checked peace (10yr)", cex.var=0.9)










#### 2. Predicted Probability & spagetti plot -- Model 17 (model with agreed partition)
data<-checked.data
sophie3<-m19
set.seed(123)

par(mfrow=c(1,1))
#m19: fortna10~partition

lwgvec <- seq(min(data$partition), max(data$partition), length.out=100)

predvec <- NULL
for(i in 1:length(lwgvec)){
  values <- c(1, lwgvec[i])
  pred <- sophie3$coefficients %*% values
  predvec <- c(predvec, pred)
}

predprob <- 1/(1+exp(-predvec))
#predprob <- -predvec



### WITH UNCERTAINTY 

draw <- mvrnorm(1000, sophie3$coefficients, vcov(sophie3))

predmat <- matrix(NA, nrow=1000, ncol=100)
for(j in 1:1000){
  predvec <- NULL
  for(i in 1:length(lwgvec)){
    values <- c(1, lwgvec[i])
    pred <- draw[j,] %*% values
    predvec <- c(predvec, pred)
  }
  predmat[j,] <- 1/(1+exp(-predvec))
}

# plotting confidence interval (spaghetti plot type)
# Box


par(mfrow=c(1,1))
plot(lwgvec, predprob, type="n", ylim=c(0,1), col=accent[5], axes=F,
     ylab="Survival of Peace", xlab="Partition", main="(M3) Predicted Probability")
box()
#X axis
axis(1, at=0, lab="Non-partition")
axis(1, at=1, lab="Partition")
#Y axis
axis(2, at=0.1, lab="Recurrence")
axis(2, at=1, lab="Peace")



for(i in 1:1000){
  points(lwgvec, predmat[i,], type="l", lwd=0.2, col=spectral[4])
}
lines(lwgvec, predprob, type="l", col=spectral[2], lwd=2)



#Rug at the top (placing rug at the top side=3, defulat;bottom)
#rug(data$partition)





#### 3. Predicted Probability & spagetti plot (model 18, agreed partition & peace)
data<-checked.data
sophie3<-m20
set.seed(123)


lwgvec <- seq(min(data$partition), max(data$partition), length.out=100)

predvec <- NULL
for(i in 1:length(lwgvec)){
  values <- c(1, lwgvec[i])
  pred <- sophie3$coefficients %*% values
  predvec <- c(predvec, pred)
}

predprob <- 1/(1+exp(-predvec))
#predprob <- -predvec



### WITH UNCERTAINTY 

draw <- mvrnorm(1000, sophie3$coefficients, vcov(sophie3))

predmat <- matrix(NA, nrow=1000, ncol=100)
for(j in 1:1000){
  predvec <- NULL
  for(i in 1:length(lwgvec)){
    values <- c(1, lwgvec[i])
    pred <- draw[j,] %*% values
    predvec <- c(predvec, pred)
  }
  predmat[j,] <- 1/(1+exp(-predvec))
}

# plotting confidence interval (spaghetti plot type)
# Box


par(mfrow=c(1,1))
plot(lwgvec, predprob, type="n", ylim=c(0,1), col=accent[5], axes=F,
     ylab="Survival of Peace", xlab="Partition", main="(M4) Predicted Probability")
box()
#X axis
axis(1, at=0, lab="Non-partition")
axis(1, at=1, lab="Partition")
#Y axis
axis(2, at=0.1, lab="Recurrence")
axis(2, at=1, lab="Peace")



for(i in 1:1000){
  points(lwgvec, predmat[i,], type="l", lwd=0.2, col=spectral[8])
}
lines(lwgvec, predprob, type="l", col=spectral[10], lwd=2)



#Rug at the top (placing rug at the top side=3, defaulat;bottom)



############### Finally, two scenarios!


##################################
##### probability of two scenarios   Model 19 (m19)
#################################


##### Generate the probability of war based on two different scenarios(partition/-partition)

data<-checked.data
sophie3<-m19
# setting up vectors to be filled
outmin <- outmin <- NULL  
set.seed(123)
draw <- mvrnorm(n=1000,sophie3$coefficients,vcov(sophie3))   
dim(draw)

# the vector of values: 1 for intercept
effectmin <- c(1,0)   
effectmax <- c(1,1)    

outmin <- draw%*%effectmin
outmax <- draw%*%effectmax

dim(outmin)
dim(outmax)

# in log scale (since it's a logit model..)
outmin<- 1/(1+exp(-outmin))
dim(outmin)
summary(outmin)

outmax <- 1/(1+exp(-outmax))
dim(outmax)
summary(outmax)

plot(density(outmin))
plot(density(outmax)) 

######### Now graphing!

#Prepare the box (outside)
plot(density(outmin), type="n", ylim=c(0,5), xlim=c(0,1), axes=F,
     ylab="Density", xlab="Probability of Peace", main="(M3)Two Scenarios of Partition and -Partition on Peace")
box()
axis(1, at=0, lab=0)
axis(1, at=0.5, lab=0.5)
axis(1, at=1, lab=1)


d<-density(outmin) #
s<-density(outmax) #

lines(d, col=spectral[11], lwd=0.8) 
polygon(d, col=spectral[9], density=30) 
abline(v=mean(outmin),b=0, col=spectral[11], lwd=0.7, lty=2)

lines(s, col=spectral[1], lwd=0.8) 
polygon(s, col=spectral[3], density=30, opacity=0.5)
abline(v=mean(outmax), b=0, col=spectral[1], lwd=0.7, lty=2)

#Now legend!
legend("topleft",c("Partition","Non-partition"),col=c(spectral[3],spectral[9]),lty=c(1,1), lwd="3")
plot(s)



##################################
##### probability of two scenarios   Model 20 (m20)
#################################


##### Generate the probability of war based on two different scenarios(partition/-partition)

data<-checked.data
sophie3<-m20
# setting up vectors to be filled
outmin <- outmin <- NULL  
set.seed(123)
draw <- mvrnorm(n=1000,sophie3$coefficients,vcov(sophie3))   
dim(draw)

# the vector of values: 1 for intercept
effectmin <- c(1,0)   
effectmax <- c(1,1)    

outmin <- draw%*%effectmin
outmax <- draw%*%effectmax

dim(outmin)
dim(outmax)

# in log scale (since it's a logit model..)
outmin<- 1/(1+exp(-outmin))
dim(outmin)
summary(outmin)

outmax <- 1/(1+exp(-outmax))
dim(outmax)
summary(outmax)

plot(density(outmin))
plot(density(outmax)) 

######### Now graphing!

#Prepare the box (outside)
plot(density(outmin), type="n", ylim=c(0,4), xlim=c(0,1), axes=F,
     ylab="Density", xlab="Probability of Peace", main="(M4) Two Scenarios of Partition and -Partition on Peace")
box()
axis(1, at=0, lab=0)
axis(1, at=0.5, lab=0.5)
axis(1, at=1, lab=1)


d<-density(outmin) #
s<-density(outmax) #

lines(d, col=spectral[11], lwd=0.8) 
polygon(d, col=spectral[9], density=30) 
abline(v=mean(outmin),b=0, col=spectral[11], lwd=0.7, lty=2)

lines(s, col=spectral[1], lwd=0.8) 
polygon(s, col=spectral[3], density=30, opacity=0.5)
abline(v=mean(outmax), b=0, col=spectral[1], lwd=0.7, lty=2)

#Now legend!
legend("topleft",c("Partition","Non-partition"),col=c(spectral[3],spectral[9]),lty=c(1,1), lwd="3")
plot(s)







#Save everything!

setwd("C://Users/Sophie/Desktop/2013 spring/PS 659 Civil War/Final_replication/images")


################ Bayesian Analysis

#create the dataset


lmodel1<-zelig(formula= survivalofpeace~partition + tag(1|ethnic), data=datawinner, model="logit.mixed") 
summary(lmodel1)

